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to treatment in terms of uniform contractions. This approach provides an interesting 
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Chapter 1 
Introduction 



The art of doing mathematics consists in finding that special 
case which contains all the germs of generality. 

- David Hilbert 

Real world complex systems can be viewed and modeled as networks of interacting 
elements HHUE]. Examples range from geology [4| and ecosystems to math- 
ematical biology H and neuroscience Q as well as physics of neutrinos [8| and 
superconductors [9 |. Here we distinguish the structure of the network, the nature of 
the interaction, and the (isolated) dynamical behavior of individual elements. 

During the last fifty years, empirical studies of real complex systems have led to 
a deep understanding of the structure of networks, interaction properties, and iso- 
lated dynamics of individual elements, but a general comprehension of the resulting 
network dynamics remains largely elusive. 

Among the large variety of dynamical phenomena observed in complex net- 
works, collective behavior is ubiquitous in real world networks and has proven to be 
essential to the functionality of such networks lfT0l[TTl[T2l[T3ll . Synchronization is 
one of the most pervasive form collective behavior in complex systems of interacting 
components lfl4l [TBI [161 [171 [181 . Along the riverbanks in some South Asian forests 
whole swarms of fireflies will light up simultaneously in a spectacular synchronous 
flashing. Human hearts beat rhythmically because thousands of cells synchronize 
their activity fl4l . while thousands of neurons in the visual cortex synchronize their 
activity in response to specific stimuli 1191 . Synchronization is rooted in human life, 
from the metabolic processes in our cells to the highest cognitive tasks G0ll2"Tl . 

Synchronization emerges from the collaboration and competition of many ele- 
ments and has important consequences to all elements and network functioning. 
Synchronization is a multi-disciplinary discipline with broad range applications. 
Currently, the field experiences a vertiginous growth and significant progress has 
already been made on various fronts. 

Strikingly, in most realistic networked systems where synchronization is rele- 
vant, strong synchronization may also be related to pathological activities such as 
epileptic seizures El and Parkinson's disease G3ll in neural networks, to extinction 
in ecology [24|, and social catastrophes in epidemic outbreaks [25 1. Of particular 
interest is how synchronization depends on various structural parameters such as 
degree distribution and spectral properties of the graph. 
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1 Introduction 



In the mid-nineties Pecora and Carroll [26] have put forward a paradigmatic 
model of diffusively coupled identical oscillators on complex networks. They have 
shown that complex networks of identical nonlinear dynamical systems can glob- 
ally synchronize despite exhibiting complicated dynamics at the level of individual 
elements. 

The analysis of synchronization in complex networks has benefited from ad- 
vances in the understanding of the structure of complex networks 112711281 . Barahona 
and Pecora ||29ll have shown that well-connected networks - with so-called small- 
world structure - are easier to globally synchronize than regular networks. Motter 
and collaborators BUI have shown that heterogeneity in the network structure hin- 
ders global synchronization. These results form only the beginning of a proper un- 
derstanding of the connections between network structure and the stability of global 
synchronization. 

The approach put forward by Pecora and Carroll, which characterized the sta- 
bility of global synchronization, is based on elements of the theory of Lyapunov 
exponents |3 lj. The characterization of stability via theory of Lyapunov exponents 
has many additional subtleties, in particular, when it comes to the persistence of 
stability under perturbations. Positive solution to the persistence problem requires 
the analysis of the so called regularity condition, which is tricky and difficult to 
establish. 

We consider the fully diffusive case - the coupling between oscillators depends 
only on their state difference. This model is amenable to full analytical treatment, 
and the stability analysis of the global synchronization is split into contributions 
coming solely from the dynamics and from the network structure. The stability con- 
ditions in this case depend only on general properties of the oscillators and can be 
obtained analytically if one possesses knowledge of global properties of the dynam- 
ics such as the boundedness of the trajectories. We establish the persistence under 
nonlinear perturbations and linear perturbations. Many conclusions guide us toward 
the ultimate goal of understanding more general collective behavior 



Chapter 2 

Graphs : Basic Definitions 



Can the existence of a mathematical entity be proved without 
defmiing it? 

- Jacques Hadamard 



2.1 Adjacency and Laplacian Matrices 

A network is a graph G comprising a set of N nodes (or vertices) connected by 
a set of M links (or edges). Graphs are the mathematical structures used to model 
pairwise relations between objects. We shall often refer to the network topology, 
which is the layout pattern of interconnections of the various elements. Topology 
can be considered as a virtual shape or structure of a network. 

The networks we consider here are simple and undirected. A network is called 
simple if the nodes do not have self connections, and undirected if there is no dis- 
tinction between the two vertices associated with each edge. A path in a graph is a 
sequence of connected (non repeated) nodes. From each of node of a path there is a 
link to the next node in the sequence. The length of a path is the number of links in 
the path. See further details in Ref. |32|. 

For example, lets consider the network in Fig. |2.1fr ). Between to the nodes 2 and 
4 we have three paths {2, 1,3,4}, {2,5,3,4} and {2,3,4}. The first two have length 
3, and the last has length 2. Therefore, the path {2,3,4} is the shortest path between 
the node 2 and 4. 

The network diameter d is the greatest length of the shortest path between any 
pair of vertices. To find the diameter of a graph, first find the shortest path between 
each pair of vertices. The greatest length of any of these paths is the diameter of the 
graph. If we have an isolated node, that is, a node without any connections, then we 
say that the diameter is infinite. A network of finite diameter is called connected. 

A connected component of an undirected graph is a subgraph with finite diameter. 
The graph is called directed if it is not undirected. If the graph is directed then there 
are two connected nodes say, u and v, such that u reachable from v, but v is not 



reachable from u. See Fig. 2.1 for an illustration. 

The network may be described in terms of its adjacency matrix A, which encodes 
the topological information, and is defined as 



A,j = 



1 if i and j are connected 
otherwise . 
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2 Graphs : Basic Definitions 




connected disconnected directed 

Fig. 2.1 Examples of undirected a) and b) and directed c) graphs. The diameter of graph a) is d = 2 
hence, the graph is connected. The graph b) is disconnected, there is no path connecting the nodes 
1 and 2 to the remaining nodes, the diameter is d = °°. However, the graph has two connected 
components, the upper (1,2) with diameter d = l, and the lower nodes (3,4,5) with diameter d = 2. 
Graph c) is directed, the arrow tells the direction of the connection, so node 1 is reachable from 
node 2, but not the other way around. 



An undirected graph has a symmetric adjacency matrix. The degree ki of the ith 
node is the number of connections it receives, clearly 

h = ^Ajj. 

J 

Another important matrix associated with the network is the combinatorial 
Laplacian matrix L, defined as 

ki if i = j 

Lij = t — 1 if « and j are connected 
otherwise . 

The Laplacian L is closely related to the adjacency matrix A. In a compact form 
it reads 

L=D -A, 



where D = diag(£i, • • • is the matrix of degrees. We depict in Fig. 2.2 distinct 
networks of size 4 and their adjacency and laplacian matrices. 



2.2 Spectral Properties of the Laplacian 

The eigenvalues and eigenvectors of A and L tell us a lot about the network struc- 
ture. The eigenvalues of L for instance are related to how well connected is the 
graph and how fast a random walk on the graph could spread. In particular, the 
smallest nonzero eigenvalue of L will determine the synchronization properties of 
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Fig. 2.2 Networks containing four nodes. Their adjacency and laplacian matrices are represented 
by A and L. Further details can be found in Table|2.2] 



the network. Since the graph is undirected the matrix L is symmetric its eigenvalues 
are real, and L has a complete set of orthonormal eigenvectors [9] The next result 
characterizes important properties of the Laplacian 

Theorem 1 Let G be an undirected network and L its associated Laplacian. Then: 

a) L has only real eigenvalues, 

b) is an eigenvalue and a corresponding eigenvector is 1 = (1 , 1 , • • • , 1 )*, where * 
stands for the transpose. 

c) L is positive semidefinite, its eigenvalues enumerated in increasing order and 
repeated according to their multiplicity satisfy 

= X\ < A 2 < • • • < K 

d) The multiplicity ofO as an eigenvalue of I, equals the number of connect compo- 
nents of G. 

Proof : The statement a) follows from the fact that L is symmetric L = L* , see 
Ap . [A| Theorem [9] To prove b) consider the / = (1, 1)* and note that 



(Ll) i = 1 £Lij = k i - 1 £A i j = 

j 



(2.1) 
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2 Graphs : Basic Definitions 



Item c) follows from the Gershgorin theorem, see Ap. [A|Theorem[8] The nontrivial 
conclusion d) is one of the main properties of the spectrum. To prove the statement 
d) we first note that if the graph G has r connected components G\ , • • • , G r , then is 
possible to represent L such that it splits into blocks L\ . ■ ■ ■ ,L r . 

Let m denote the multiplicity of 0. Then Each L, has an eigenvector z, with 
as an eigenvalue. Note that z; = (zj , ■ ■ -z") can be defined as zj is equal to 1 if j 
belongs to the component i and zero otherwise, hence m > r. It remains to show that 
any eigenvector g associated with is also constant. Assume thatg is a non constant 
eigenvector associated with 0, and let g( > be the largest entry of g. Then 

( L g)e = Y, L ijgj 

j 

j 

since g is associated with the eigenvalue we have 

8i = , ■ 

This means that the value of the component g( is equal to the average of the values 
assigned to its neighbors. Hence g must be constant, which completes the proof. □ 
Therefore, X2 is bounded away from zero whenever the network is connected. 
The smallest non-zero eigenvalue is known as algebraic connectivity, and it is often 
called the Fiedler value. The spectrum of the Laplacian is also related to some other 
topological invariants. One of the most interesting connections is its relation to the 
diameter, size and degrees. 

Theorem 2 Let G be a simple network of size n and L its associated Laplacian. 
Then: 

1. fi33$h> 

nd 

2. fiEHX 2 < -^-rki 

n — 1 

We will not present the proof of the Theorem here, however, they can be found in 
references we provide in the theorem. We suggest the reader to see further bounds 
on the spectrum of the Laplacian in Ref. ||35l . Also Ref. (36 1 presents many appli- 
cations of the Laplacian eigenvalues to diverse problems. One of the main goals in 
spectral graph theory is the obtain better bounds by having access to further infor- 
mation on the graphs. 

For a fixed network size, the magnitude of A2 reflects how well connected is 
graph. Although the bounds given by Theorem [2] are general they can be tight for 
certain graphs. For the ring the low bound on A2 is tight. This implies that as the size 
increase - consequently also its diameter - X2 converges to zero, and the network 
becomes effectively disconnected. In sharp contrast we find the star network. In this 
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case, the upper bound in i) is tight. The star diameter equals to two, regardless the 
size and X2 — 1. See the table for the precise values. 



Table 2.1 Network of n nodes. Examples of such networks are depicted in Fig. 



2.2 



Network X2 k„ k\ D 

Complete n re— In— 1 1 

, . (re + l)/2ifreisodd 

ring 2-2cos — 2 2 v , " 

\n J re/2 11 n is even 

Star 1 n - 1 1 2 



The networks we encounter in real applications have a wilder connection struc- 
ture. Typical examples are cortical networks, the Internet, power grids and metabolic 
networks [ 1|. These networks don't have a regular structure of connections such as 
the ones presented in Fig. 2.2 We say that the network is complex if it does not 
possess a regular connectivity structure. 

One of the goals is the understand the relation between the topological organi- 
zation of the network and its relation functioning such as its collective motion. In 
Fig. 2.3 we depict two networks used to model real networks, namely the Barabasi- 
Albert and the Erdos-Renyi Networks. 




Scale-Free Erdos-Renyi 

Fig. 2.3 Some examples of complex networks. 

The Erdos-Renyi network is generated by setting an edge between each pair of 
nodes with equal probability p, independently of the other edges. If p ^> \nn/n, 
then a the network is almost surely connected, that is, as tends to infinity, the 
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2 Graphs : Basic Definitions 



probability that a graph on n vertices is connected tends to 1. The degree is pretty 
homogeneous, almost surely every node has the same expected degree 11271 . 

The Barabasi-Albert network possesses a great deal of heterogeneity in the 
node's degree, while most nodes have only a few connections, some nodes, termed 
hubs, have many connections. These networks do not arise by chance alone. The 
network is generated by means of the cumulative advantage principle - the rich gets 
richer. According to this process, a node with many links will have a higher prob- 
ability to establish new connections than a regular node. The number of nodes of 
degree k is proportional to k~P. These networks are called scale-free networks (l). 
Many graphs arising in various real world networks display similar structure as the 
Barabasi-Albert network EJEl- 



Chapter 3 

Nonlinear Dynamics 



How can intuition deceive us at this point ? 

- Henri Poincare 



Let D be an open simply connected subset of M. m , m > 1, and let / £ C(D,W") for 
some r > 2. We assume that the differential equation 

(3-D 

models the dynamics of a given system of interest. Now since / is differentiable 
the Picard-Lindelof Theorem guarantees the existence of local solutions, see Ap. 



[B] Theorem 16 We wish to guarantee that the solutions also exist globally. This 
requires further hypothesis on the behavior of the vector field. We are interested in 
systems that dissipate the volumes of M m - called dissipative systems. 



3.1 Dissipative Systems 



We say that set Q C M." 1 under the dynamics of Eq. ( 3. 1 1 is positively invariant if the 
trajectories starting at the set never leave it in the future, that is, if jc(fo) £ Q then 
x(t) £ £2 for all t > to Intuitively, it means that once the trajectory enters Q it never 
leaves it again. The system is called dissipative if the solutions enter a positively 
invariant set £2 C D in finite time. Q is called absorbing domain of the system. 
The existence of an absorbing domain guarantees that the solutions are bounded, 



hence, the extension result in Ap. [BjTheorem 17 assures the global existence of the 
solutions. 

The question is then how to obtain the absorbing domains. Note that whenever / 



is nonlinear finding the solutions of Eq. 3.1 can be a rather intricate problem. And 



usually we won't be able to do it analytically. So we need new machinery to address 
the problem on absorbing domains. A method by Lyapunov allows us to obtain such 
domains without finding the trajectories. The technique infers the existence of the 
absorbing domains in relation to some properties of a scalar function - the Lyapunov 
function. 
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3 Nonlinear Dynamics 



We will study notions relative to connected nonempty subsets Q of R m . A func- 
tion V : W" —> M is said to be positive definite with respect to the set B if V (x) > 
for all x G M.^Q. It is radially unbounded if 

lim V(x) = oo. 

\\x\\— 

Note that this condition guarantees that all level sets of V are bounded. This fact 
plays a central role in the analysis. We also define V' : K m — > K as 

V'(x)=VV(jc) •/(*). 

where • denotes the Euclidean inner product. This definition agrees with the time 
derivative along the trajectories. That is, if x(t ) is a solution of Eq. ( 3.1 1, then by the 
chain rule we have 

dV{x(t )) 
dt 

the main result is then the following 



= V'(x(t)). 



Theorem 3 (Lyapunov) Let V : W" — > E be radially unbounded and positive defi- 
nite with respect to the set Q. C D. Assume that 

V'(x) <0forallxeR'"\£2 



Then all trajectories of Eq. (3.1 \ eventually enter the set Q, in other words, the 
system is dissipative. 

Proof: Note that for any trajectory x(t) in virtue of the fundamental theorem of 
the calculus 



V(x(t))-V(x(s)) = J V'{x{u))du 



<0. 



So V(x(t)) < V(x(s)) for any t > s, and V is decreasing along solutions and is 
radially unbounded, the level sets 

S a = {xeR m :V(x)<a} 

are positively invariant. Hence, the solutions are bounded, and will lie in smaller 
level sets as time increase until the trajectory enters Q. It remains to show that once 
the solutions lie in £2 they don't leave it. 

Suppose x(t) leaves Q at to and let b = V(x(to)). The level set Sb is closed, and 
there is a ball B r (x(to)) such that x(to + e) G B r (x(to))\Sb for some small e. Hence, 
V(x(to + e)) > V (jc(fo)) contradicting the fact that V is decreasing along solutions. 
□ 

There are also converse Lyapunov theorems [37|. Typically if the system is dis- 
sipative (and have nice properties) then there exists a Lyapunov function. Although 
the above theorem is very useful, since we don't need knowledge of the trajectories, 
the drawback is the function V itself. There is no recipe to obtain a function V ful- 
filling all these properties. One could always try to guess the function, or go for a 



3.2 Chaotic Systems 1 1 

general form such as choosing a quadratic function V. We assume that the Lyapunov 
function is given. 

Assumption 1 There exists a symmetric positive matrix Q such that 

V(x) = ^(x-afQ(x-a). 
where a € W". Consider the set Q := {x e K m | (x -o) t Q(x -a) < p 2 }, then 

V'(x) <o,Vxem. m \£2. 

Under Assumption [T] Theorem [5] guarantees that £2 is positively invariant and 
that the trajectories of Eq. ( |3.1| l eventually enter it. So, Q is the absorbing domain 
of the problem. The solutions are, therefore, globally defined. 



3.2 Chaotic Systems 

Since the system Eq. ( |3.1| l is dissipative the solutions accumulates in a neighbor- 
hood of abounded set A C Q. The set A is called attractor. We focus on the situation 
where A is a chaotic attractor. Now, the definition of a chaotic attractor is rather in- 
tricate - there is even a general definition, the important properties for us are that so- 
lutions on the attractor are aperiodic, i.e., there is no T > such that x(t) =x(t + 1), 
and the solutions exhibits sensitive dependence on initial conditions. Sensitive de- 
pendence on initial conditions means that nearby trajectories separate exponentially 
fast. 




Fig. 3.1 If jc(0) ^ y(0) are both are in a neighborhood of the attractor, then for small times \\x(t ) — 
y(t)\\ « ||x(0) -.y(0)j|e c '',for some b > 0. 



If the system is chaotic no matter how close two solutions start, they move apart 
when they are close to the attractor. Hence, arbitrarily small modifications of ini- 
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3 Nonlinear Dynamics 



tial conditions typically lead to quite different states for large times. This sensitive 
dependence on initial conditions is one of the main features of a chaotic system. 
Exponential divergency cannot go on for ever, since the attractor is bounded, it is 
possible to show that the trajectories will come close together in the future [ 38 ]. 



3.2.1 Lorenz Model 

The Lorenz model exhibits a chaotic dynamics (39). Using the notation 




the Lorentz vector field reads 




/(*) 



where we choose the classical parameter values a = 10, r = 28, b = 8/3. For these 
parameters the Lorenz system fulfills our assumption [T]on dissipativity. 

Proposition 1 The trajectories of the Lorenz eventually enter the absorbing domain 

h 2 r 2 

3 . „JZ i -.,,2 i „r„ o„\2 ° r 



Q = <x e K J : rx z + oy z + o{z - 2r) z < 



b-l 

Proof: Consider the function 

V(x) = {x-a)^Q{x-a) 

where a = (0,0, 2r) and Q = diag(r, (7,(7), note that the matrix is positive-definite. 
The goal is to find a bounded region Q - defined by means of a level set of V - such 
that V' < in the exterior of Q and then apply Theorem[3] To this end, we compute 
the derivative, 

V'{x) =2(x-aYQf(x) 

= 2orx(y — x) + Gyx(r-z) — (7y 2 — baz{z — 2r) + o(z — 2r)xy 
= —2orx 2 — ay 2 + oyx(r — z) — boz(z — 2r) — a(r — z)xy 
= -2a[rx 2 +y 2 +b{z-r) 2 -br 2 ]. 

Consider the ellipsoid E defined by rx 2 +y 2 + b(z — r) 2 < br 2 , hence, in the 
exterior of E we have V'. Now we take c to be the largest value of V in E, and we 
define Q = {x € 1R 3 : V(x) < c}. The solutions will eventually enter Q and remain 
inside since V' < in the exterior of Q, and once the trajectory enters in Q it never 
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leaves the set. It remains to obtain the parameter c. This can be done by means of a 
Lagrande multiplier. After a computation - see Appendix C of Ref. [40 1 - we obtain 
c = b 2 r z /(b-l)forb>2, and a > 1. □ 

Inside the absorbing set Q the trajectory accumulates on the chaotic attractor. 
We have numerically integrated the Lorentz equations using a fourth order Runge- 
Kutta, the initial conditions are x(0) = —10, y(0) = 10, z(0) = 25. We observe that 
the trajectory accumulates on the so called Butterfly chaotic attractor ll39l . see Fig. 

EH 



Fig. 3.2 The trajectories of the Lorenz system eventually enters an absorbing domain and accumu- 
lates on a chaotic attractor. This projection of attractor resembles a butterfly - the common name 
of the Lorenz attractor. 



Close to the attractor nearby trajectories diverge. To see this phenomenon in 
a simulation let us consider a distinct initial condition JE(0) = (Jc(0),y(0),z(0))*. 
We consider f(0) = -10.01, y(0) = 10, z(0) = 25. Note that the initial difference 
||jt(0) — jc (0) 1 1 2 = 0.01 becomes as large as the attractor size in a matter of 6 cycles, 
see Fig. |3.3| 



3.3 Diffusively Coupled Oscillators 

We introduce now the network model. On the top of each node of the network we 



introduce a copy of the system Eq. (3.1 1. Then the influence that the neighbor j 



exerts on the dynamics of the node i will be proportional to the difference of their 
state vector Xj(t) —Xi(t), This type of coupling is called diffusive - it tries to equate 
to state of the nodes. 

We label the nodes according to their degrees k„ > ■ ■ ■ > k^ > k\, where k\ and k n 
denote the minimal and maximal degree, respectively. The dynamics of a network 
of n identically diffusively coupled elements is described by 
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3 Nonlinear Dynamics 




Fig. 3.3 Two distinct simulations of the time series x(t) and x(t) of the Lorentz systems. The 
difference between the trajectories is of 0.01, however this small difference grows with time until 
a point where the difference is as large as the attractor itself. 



dXi 
dt 



--f(Xi) + aY,Aij(Xj-Xi), 

7=1 



(3.2) 



where a the overall coupling strength. In Eq. ( 3.2 1 the coupling is given in terms of 
the adjacency matrix. We can also represent the coupled equations in terms of the 
network laplacian. Consider the coupling term 



^•^■17 (■*•/' •*•') — 7 . -Aj jX j kjX, 
7=1 y'=i 
n 

= E ( A U - 8 iJ k ') x j 
7=1 

where is the Kronecker delta, and recalling that L,- ; = —A[j we obtain Hence, 
the equations read 



dXi 
dt 



--f(xi) -a Y, L iJ x J- 



(3.3) 



7=1 



The dynamics of such diffusive model can be intricate. Indeed, even if the iso- 
lated dynamics possesses a globally stable fixed point the diffusive coupling can 
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lead to the instability of the fixed points and the systems can exhibit an oscillatory 
behavior. Please, see |41| for a discussion and further material. We will not focus 
on such scenario of instability, but rather on how the diffusive coupling can lead to 
synchronization. 

Note that due to the diffusively nature of the coupling, if all oscillators start with 
the same initial condition the coupling term vanishes identically. This ensures that 
the globally synchronized state 

Xl(t) =xz(t) =---=x„(t)=s(t), 

is an invariant state for all coupling strength a. The question is then the stability of 
synchronized solutions, which takes place due to coupling. Note that, if a — the 
oscillators are decoupled, and Eq. ( |3.3| > describes n copies of the same oscillator with 
distinct initial conditions. Since, the chaotic behavior leads to a divergence of nearby 
trajectories, without coupling, any small perturbation on the globally synchronized 
motion will grow exponentially fast, and lead to distinct behavior between the node 
dynamics. 

The correct way to see the invariant of globally synchronized motion is as fol- 
lows. First consider 

X = co\(x u --- ,X n ), 

where col denotes the vectorization formed by stacking the columns vectors into 
a single column vector. Similarly 

F(X)=col(f(xi),-J(x n )), 

then Eq. ( |3.3| l can be arranged into a compact form 

d £=F{X)-a(L®I m )X (3.4) 

where ® is the Rronecker product, see Appendix [A| Let ( P(-,t) be the flow of 
Eq. (3.4 1, the solution of the equation with initial condition Xq is given by X(t) = 
<P(XqJ). Consider the synchronization manifold 

.£ = {xi e R m : Xi{t) =s(t) for 1 < i < «}, 

then we have the following result 

Proposition 2 is an invariant manifold under the flow ) 

Proof: Recall that 1 g R" is such that every component is equal to 1. Let X(t) — 
l®s(t), note that 

dX(t) ds(t) 
dt dt 
We claim that X(t) is a solution of the equations of motion. 
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3 Nonlinear Dynamics 



dX ^=F(X{t))-a{L®I m )X{t) 
at 

= F{l®s{t))-a{L®l m )l®s{t) 
where in the last passage we used Theorem () together with L 1 = and F(l i 



s(t)) = ltg>f(s(t)). By the Picard-Lindelof Theorem 16 we have that that X(t) = 
<P(1 ®s(Q),t) € 'Jl for all t. □ 

If the oscillators have the same initial condition, their evolution will be exactly 
the same forward in time no matter the value of the coupling strength. 

In the above result we have looked at the network not as coupled equation but as 
a single system in the full state space M. mn . We prefer to keep the picture of coupled 
oscillators. These pictures are equivalent and we interchange them whenever it suits 
our purposes. The important questions are 

Boundedness of the solutions Xi(t). 

Stability of the globally synchronized state (synchronization manifold). 

We wish to address the local stability of the globally synchronized state. That is, 
if all trajectories start close together ||jc,-(0) — acy(0)|| < £, for any i and j and some 
small e, would they converge to in other words, would 

lim \\xj(t) —Xj(t)\\ = 



or would the trajectories split apart? The goal of the remaining exposition is to pro- 
vide positive answers to these questions. To this end, we review some fundamental 
results needed to address such points. 



Chapter 4 

Linear Differential Equations 



The more you know, the less sure you are. 

- Voltaire 



The question concerning the local stability of a given trajectory s(t) leads to the 
stability analysis of the trivial solution of a nonautonomous linear differential equa- 
tion. The analysis of the dynamics in a neighborhood of the solutions is performed 
by using the variational equation. The trajectory s(t) is stable when the stability of 
the trivial solution of the variational equation is preserved under small perturbations. 



4.1 First Variational Equation 

Letj(O) be close to s(0). Each of these distinct points has its behavior determined 
by the equation of motion Eq. ([XT}. We can follow the dynamics of the difference 

z(t)=y(t)-s(t) 

which leads to the variational equations governing its evolution 

dz{t) 



dt 



=m))-M*)) 



= f(s(t)+z(t))-f(s(t)), 

now since ||z|| is sufficiently small we may expand the function /in Taylor series 

f( S (t) +z(t)) =f(s(t)) +Df(s(t))z(t) +R(z(t)) 

where Df(s(t)) along the trajectory s(t), and by the Lagrange theorem ll42ll 

||/?(z(f))|| = 0(||z(f)|| 2 ). 

Truncating the evolution equation of z, up to the first order, we obtain the first vari- 
ational equation 
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4 Linear Differential Equations 



Note that the above equation is non-autonomous and linear. Moreover, since s(t) 
lies in a compact set and / is continuously differentiable, by Weierstrass Theorem 
11421 . Df(s(t)) is a bounded matrix function. If ||z(f)|| — > the two distinct solutions 
converge to each other and have an identical evolution. 

The first variational equation plays a fundamental role to tackling the local sta- 
bility problem. Suppose that somehow we have succeeded to demonstrate that the 
trivial solution of the first variational equation is stable. Note that this does not com- 
pletely solve our problem, because the Taylor remainder acts as a perturbation of the 
trivial solution. Hence, to guarantee that the problem can be solved in terms of the 
variational equation we must also obtain conditions on the persistence of the stabil- 
ity of trivial solution under small perturbation. There is a beautiful and simple, yet 
general, criterion based on uniform contractions. We follow closely the exposition 
inRef. EHUD. 



4.2 Stability of Trivial Solutions 

Consider the linear differential equation 

dx , . 

#=u<fyx (4.i) 

where U(t) is a continuous bounded linear operator on W for each t > 0. 



The point* = is an equilibrium point of the equation Eq. (4.1 1. Loosely speak- 
ing, we say an equilibrium point is locally stable if the initial conditions are in a 
neighborhood of zero solution remain close to it for all time. The zero solution is 
said to be locally asymptotically stable if it is locally stable and, furthermore, all 
solutions starting near tend towards it as t 



The time dependence in Eq. (4.1 1 introduces of additional subtleties B31 . There- 



fore, we want to state some precise definitions of stability 

Definition 1 (Stability in the sense of Lyapunov) The equilibrium point x* = is 
stable in the sense of Lyapunov at t — to if for any £ > there exists a <5(fo,e) > 
such that 

||x()o)|| <S=H|x(f)|| <£, Vf >t 

Lyapunov stability is a very mild requirement on equilibrium points. In particular, 
it does not require that trajectories starting close to the origin tend to the origin 
asymptotically. Also, stability is defined at a time instant tQ. Uniform stability is a 
concept which guarantees that the equilibrium point is not losing stability. We insist 
that for a uniformly stable equilibrium point x*, 8 in the Definition 4.1 not be a 
function of ?o, so that equation may hold for all tQ. Asymptotic stability is made 
precise in the following definition: 
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Definition 2 (Asymptotic stability) An equilibrium point x* = is asymptotically 
stable att — to if 

1. x* = is stable, and 

2. x* = is locally attractive; i.e., there exists 8 {to) such that 

||x(f )|| < 8 =*> limx(f) =0 

Definition 3 (Uniform asymptotic stability) An equilibrium point x* = is uni- 
form asymptotic stability if 

1. x* = is asymptotically stable, and 

2. there exists 8q independent of to for which equation holds. Further, it is required 
that the convergence is uniform. That is, for each £ > a corresponding T = 
T{e) >0such thatif\\x{s)\\ < 8o for some s>0 then \\x{t)\\ < £ for all t >s + T. 

We shall focus on the concept of uniform asymptotic stability. To this end, we 
wish to express the solutions of the linear equation in a closed form. The theory of 
differential equations guarantees that the unique solution of the above equation can 
be written in the form 

x{t)=T{t,s)x{s) 

where T{t,s) is the associated evolution operator [44 1. The evolution operator satis- 
fies the following properties 

T{t 7 s)T{s,u) = T{t,u) 
T{t,s)T{s,t)=I m . 

The following concept plays a major role in these lectures 



Definition 4 Let T{t,s) be the evolution operator associated with Eq. (4.1 1. T{t,s) 
is said to be a uniform contraction if 

||T(M)|| <Ke-^- s) . 
where K and r\ are positive constants. 

Some examples of evolution operators and uniform contractions are 



Example 1 If U is a constant matrix, then Eq. \4.1\ is autonomous, and the funda- 
mental matrix reads 

T(f, 5 )=e^' u , 

T(f,s) has a uniform contraction if, and only if all its eigenvalues have negative real 
part. 

Example 2 Consider the scalar differential equation 

x = {sinlog(f + 1) +coslog(f + 1) —b}x, 
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the evolution operator reads 

T(t,s) =exp{-b(t-s) + (t + l)sinlog(t+l) - (s + 1) sin log (5 + 1)} 

Then following holds for the equilibrium point x = 

i) Ifb<\, the equilibrium is unstable. 

ii) Ifb — 1, the equilibrium is stable but not uniformly stable. 

Hi) If\<b< \Jl, the equilibrium is asymptotically stable but not uniformly stable 
or uniformly asymptotically stable. 

iv) Ifb = v2. the equilibrium is asymptotically stable. Though it is uniformly stable, 
it is not uniformly asymptotically stable. 

v) Ifb > y2, the equilibrium is uniformly asymptotically stable. 



We will show that the trivial solution of Eq. 4. 1 is uniformly asymptotically stable 
if, and only if, the evolution operator is a uniform contraction, that is, the solutions 
converge converges exponentially fast to zero. 



Theorem 4 The trivial solution of Eq. (4.1 \ is uniformly asymptotic stable if and 
only if the evolution operator is a uniform contraction. 

Proof: First suppose the evolution operator is a uniform contraction then 

||x(f)|| = \\T(t,s)x(s)\\ 

<\\T(t,s)\\\\x(s)\\ 
||x(a)||. 

Now let £ > be given, clearly if t > T, where T = T(e) is large enough then the 
\\x(t)\\ < e. Let ||*(*)|| < 5, we obtain ||jc(f)|| < Ke'^'^S <e, which implies that 

T = Tie) = —In — , 
a e 

completing the first part. 

To prove the converse, we assume that the trivial solution is uniformly asymp- 
totically stable. Then there is 8 such that for any e and T = T(e) such that for any 
IkMII < 5 we have 

\\x(t)\\<e, 

for any t > s + T. Now take E = S/k, and consider the sequence t„=s + nT. 
Note that 

||r(M)*MII<p 

for any ||jc(s)||/5 < 1, we have the following bound for the norm 

||r(M)|| = sup ||r(r,*)«|| < \. 

ii«ii<i K 
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Remember that T(t,u)T(u,s) = T(t,s). Hence, 

||r(fe,i)|| - \\T(s + 2T,s + T)T{s + T,s)\\ 

< \\t(s+2t,s+t)\\\\t(s+t,s)\\ 
1 



Likewise, by induction 
take a = Ink/T, therefore, 



\T(tn,s)\\ < 



\\T(t n ,s)\\<e- a ^- s \ 



Consider the general case t = s + u + nT, where < u < T, then the same bound 
holds 

\\T(t,s)\\<e- nTa 

< Ke-^ a , 

where K < e aT , and we conclude the desired result. □ 



4.3 Uniform Contractions and Their Persistence 

The uniform contractions have a rather important roughness property, they are not 
destroyed under perturbations of the linear equations. 

Proposition 3 Suppose U(?) is a continuous matrix function on M + and consider 



Eq. {4.1 1. Assume the fundamental matrix T(t,s) has a uniform contraction. Con- 
ion V(f ) satisfying 

sup||V(t)|| = 5< 



side r a continuous matrix function V(f) satisfying 



r>0 K 
then the evolution operator T(t,s) of the perturbed equation 

g = [U(f)+V(f)]y, 
also has a uniform contraction satisfying 

||t(M)|| <*«-*'-'>, 

where y = r\ — 8K. 



22 4 Linear Differential Equations 

Proof: Let us start by noting that the evolution operator T(t,s) also satisfies the 
differential equation of the unperturbed problem 

jT(t, S ) = U(t)T(t,s), 

The evolution operator T can be obtain by the variation of parameter, see Ap. [b] 
Theorem[l8] So, 

f(t,s)=T(t,s) + J T(t,u) V(u)t(u,s)du, 
using the induce norm, for t > s, 

\\f(t,s)\\ <Ke-rt-* +8K J\-^'-^\\t{u,s)\\du. 
Let us introduce the scalar function w(u) — e -7 ^' - ") \f(t ,s)\, then 
w{t) <Kw{s)+K8 [ w{u)du, 



for all t > s. Now we can use the Gronwall's inequality to estimate w(t), see Ap.[B| 
Theorem|2] this implies 

w(t) < Kw(s)e SK ^- s \ 

consequently 

||fM|| <Ke^- KS ^- s \ 

□ 

The roughness property of uniform contraction does the job and guarantees that 
the stability of the trivial solution is maintained. The question now turns to how 
to obtain a criterion for uniform contractions. There are various criteria, and the 
following suits our purposes 



Lemma 1 (Principle of Linearization) Assume the the fundamental matrix T(t,s) 

i. Consider the 

U(0y+R(y), 



ofEq. (4.1 \ has a uniform contraction. Consider the perturbed equation 

dy 



dt 

and assume that 

||R(y)||<M||y|| 1+c , 

for some c > 0. Then the origin is exponentially asymptotically stable. 

Proof: Note that we can write ||/?(y)|| < K\\y\\ , where K = M\\y\\. Now given a 
neighborhood of the trivial solution \\y\\ < 5 is possible to control K < e. Applying 
the previous Proposition [3] we conclude the result. □ 



4.4 Criterion for Uniform Contraction 
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This result can be used to prove that if the origin of a nonlinear system is uni- 
formly asymptotically stable then the linearized system about the origin describes 
the behavior of the nonlinear system. 



4.4 Criterion for Uniform Contraction 

The question now concerns the criteria to obtain a uniform contraction. There are 
many results in this direction, we suggest Ref. Il43ll . We present a criterion that best 
suits our purpose. The criterion provides a condition only in terms of the equation, 
and requires no knowledge of the solutions. 

Theorem 5 Let U(f ) = [Uij(t)] be a bounded, continuous matrix function on W" on 
the half-line and suppose there exists a constant 7] > such that 



for all t > and i — 1 , ■ ■ ■ , m. Then the evolution operator is a uniform contraction. 

Proof: We use the norm || • ||oo and its induced norm, see ExHin Ap. [A] Let 
x(t) be a solution. For a fixed time u > and let ||jc(h)||^, = Xi{uf. Note x(t) is 
a differentiable function and the norm a continuous function Xi(t) 2 will also be the 
norm in an open interval / = (u — a,u + a) for some a > 0. Therefore, 



m 



0s(O + £|tfy(*)l<-»7<o, 



(4.2) 



7=1, 

m 



1 d 




1 d 




2dt 



2 dt 




111 



= Uu(t)xf(t)+Y J U i jX i (t)xj(t) 



7=1, 



m 



<U a (t)x}(t)+Y,\U t} (t)\x}(t), 



7=1, 



and consequently, 



2dt 



1 d 




Using the condition 
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7=1, 



replacing in the inequality 



an integration yields 



1 d 

2dt 



\x{t)\\i<-nHmt 



(4.3) 



||*(0l£<ll*(*)l£-2fi / Wxmidx 



for all t,s € I and t > s. Applying the Gronwall inequality we have which implies 

||*(0H-<«-' ?(t -' ) ||*WII- C 4 - 4 ) 

Next note that the argument does not depend on the particular component i, because 
we assume that Eq. (4.3 i is satisfied for any 1 < i < m. So the norm will satisfy the 
bound in Eq. 4.4 for any compact set of R+. Noting that all norms are equivalent in 
finite dimensional spaces the result follows □. 



Chapter 5 

Stability of Synchronized Solutions 



Things which have nothing in common cannot he understood, 
the one by means of the other; the conception of one does not 
involve the conception of the other 

— Spinoza 

We come back to the two fundamental questions concerning the boundedness of the 
solutions and the stability of the globally synchronized in networks of diffusively 
coupled oscillators. 



5.1 Global Existence of the solutions 

The remarkable property of the networks of diffusively coupled dissipative oscilla- 
tors is that the solutions are always bounded, regardless the coupling strength and 
network structure. The two main ingredients for such boundedness of solutions are: 

- Dissipation of the isolated dynamics given in terms of the Lyapunov function 
V. 

- Diffusive coupling given in terms of the laplacian matrix 

Under these two conditions we can construct a Lyapunov function for the whole 
system. The result is then the following 

Theorem 6 Consider the diffusively coupled network model 

n 

Xi=f{xi)-aY, L ijXj, 

and assume that the isolated system has a Lyapunov function satisfying Assumption 
[7] Then, for any network the solutions of the coupled equations eventually enter an 
absorbing domain Q. The absorbing set is independent of the network. 

Proof: The idea is to construct a Lyapunov function for the coupled oscillators 
in terms of the Lyapunov function of the isolated oscillators. Consider the function 
W : W m -> R where 

W(X) = l -{X-A)*{I n ®Q){X-A) 
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where X is given by the vectorization of (jci, • • - ,x„) and likewise A = (1 ® a)*, 
where again 1 = (!,■■■ ,1). The derivative of the function W along the solutions 
reads 

dW(X) 

= (X-A)*{I n ®Q)[F(X)-a{L®I m )X\ 

= (X-A)*(I n ®Q)F(X) - aX* (L®Q)X+aA* (L<g>Q) X, 



however, using the properties of the Kronecker product, see Theorem 10 and Theo- 
rem [TT| we have 

A*{L®Q) = {l®a)*{L®Q) (5.1) 
= l*L®a*Q (5.2) 

but since 1 is an eigenvector with eigenvalue we have 1*L = 0*, and consequently 

A* {L®Q)X = Q. (5.3) 

Now L is positive semi-definite and Q is positive definite, hence it follows that 
L®Q is positive semi-definite, see Theorem 14 and 

X*(L®Q)X>0. 

We have the following upper bound 

d -^<(X-Ay(I„®Q)F(X) 

= J £(x i -a)*Qf(xi) (5-4) 

i=i 

= !>'(*,•) (5.5) 

i=i 

but by hypothesis (xj — a)*Qf(xi) is negative on D\£2, hence, dW jdt is negative on 
D"\£2", since Q depends only on the isolated dynamics the result follows. □ 
This means that the trajectory of each oscillators is bounded 

ll*i(0ll<* 

where K is a constant and can be chosen to be independent of the node i and of the 
network parameters such as degree and size. 



5.2 Trivial example: Autonomous linear equations 
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5.2 Trivial example: Autonomous linear equations 

Before we study the stability of the synchronized motion in networks of nonlinear 
equations, we address the stability problem between two mutually coupled linear 
equations. The following example is pedagogic and bears all the ideas of the prove 
of the general case. Consider the scalar equation 

dx 

— = ax 
dt 

where a > 0. The evolution operator reads 

T{t : s)=e<'- S \ 

so solutions starting at xq are given by x(t) = e at XQ. The dynamics is rather simple, 
for all initial conditions xq ^ diverge exponentially fast with rate of divergency 
given by a. Consider two of such equations diffusively coupled 

ax\ + a(x2 —Xi) 
0x2 + (X{x\ —xj) 

The pain in the neck is that the solutions of the isolated system are not bounded. 
Since the equation is linear the nontrivial solution are not bounded. On the other 
hand, because the linearity we don't need the boundedness of solutions to address 
synchronization. If a is large enough the two systems will synchronize 

lim |jti(f) — JC2(*)| = 0- 

Let us introduce 

X = 

The adjacency matrix and Laplacian are 

A= (io) andL= (-i Y) 

The coupled equations can be represented as 

= [ah - a,L] X 



dx\ 

dt 
dx2 
lit 




gi 



According to the example [T] the solution reads 

X{t) = e^-^Xo. 



(5.6) 
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We can compute the eigenvalues and eigenvectors of the Laplacian L. An easy 
computation shows that 1 = (1,1)*/V2 is an eigenvector of associated with the 
eigenvalue 0, and v 2 = ( 1 , — 1 ) / V2 is an eigenvector associated with the eigenvalue 
A2 = 2. Note that with respect to the Euclidean inner product the set {l,v 2 } is an 
orthonormal basis of R 2 . 

To solve Eq. ( |5.6| l we note that if for a given matrix B we have that u is an eigen- 
vector associated with the eigenvalue X . Then the matrix C = B — aI has eigenvector 
u associated with the eigenvalue A — a. 

We can write 

X = c\l + c 2 v 2 , 

recalling that e B 'u = e^'u, and if B and C commute then e B+c = e B e c . Hence, the 
solution of the vector equation X(t ) reads 

X(t)=e^-- aL ^(c l l+c 2 v 2 ) (5.7) 
= c l e a, l+c 2 e {a - a ^ ) 'v 2 . (5.8) 

To achieve synchronization the dynamics along the transversal mode v 2 must be 
damped out, that is, lim,^ > ooC 2 e( a ~ a ^'v 2 = 0. This implies that 

aX 2 > a =*> a > ^~ 
X 2 

Hence, the coupling strength has to be larger than the rate of divergence of the 
trajectories over the spectral gap. This is a general principle in diffusively networks. 



5.3 Two coupled nonlinear equations 



Let us consider now the stability of two oscillators diffusively coupled. At this 
time we perform the stability analysis without using the Laplacian properties. This 
allows a simple analysis and provides the condition for synchronization in the same 
spirit as we shall use later on. 



We assume that the nodes are described by Eq. ( 3. 1 1. In the simplest case of two 



diffusively coupled in all variables systems the dynamics is described by 

--f(xi) + a(x 2 -xi) 
--f(x 2 ) + a(xi -x 2 ) 
where a is the coupling parameter. Again, note that 



dx\ 
dt 

dx 2 
dt 



x\{t) =x 2 (t) 



5.3 Two coupled nonlinear equations 
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defines the synchronization manifold and is an invariant subspace of the equations of 
motion for all values of the coupling strength. Note that in the subspace the coupling 
term vanishes, and the dynamics is the same as if the systems were uncoupled. 
Hence, we do not control the motion on the synchronization manifold. If the isolated 
oscillators possess a chaotic dynamics, then the synchronized motion will also be 
chaotic. 

Again, the problem is then to determine the stability of such subspace in terms 
of the coupling parameter, the coupling strength. It turns out that the subspace it is 
stable if the coupling is strong enough. That is, the two oscillators will synchronize. 
Note that when they synchronize they will preserve the chaotic behavior. 

To determine the stability of the synchronization manifold, we analyze the dy- 
namics of the difference z=x\ — x 2 - Our goal is to obtain conditions such that 

limz = 0, 

hence, we aim at obtaining the first variational for z. 

dz{t) _ dx\(t) _ dx 2 {t) ^ 9 ^ 

dt dt dt 

= f( Xl )-f(x 2 )-2az (5.10) 

Now if ||z(0)|| -C 1, we can obtain the first variational equation governing the per- 
turbations 

d ^± = [Df(x 1 (t))-2aI]z. (5.11) 

The solutions of the variational equation can be written in terms of the evolution 
operator 

z(t)=T{t,s)z(s) 

Applying Theorem[5]we obtain conditions for the evolution operator to possesses 
a uniform contraction. Let us denote the matrix Df(x\ (f)) = [Df(x\ (0)y]"/=i' Uni- 
form contraction requires 

m 

Df( Xl (t)) tt -2a+ £ \Df(xi(t))ij\<0 (5.12) 



for all t > 0, similarly 



a c = sup 

xe£i. 

\<i<m 



£ |zy(*(0)y|+iy(x(0)s 



(5.13) 



since Q is limited and connected in virtue of the Weierstrass Theorem a c exists. 
Note that a c is closely related to the norm of the Jacobian ||D/(jt)||oo. Interestingly, 
a c can be computed only by accessing the absorbing domain and the Jacobian. Note 
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that this bound for critical coupling is usually larger than needed to observe syn- 
chronization. However, this bound is general and independent of the trajectories, 
and guarantee a stable and robust synchronized motion. 

The trivial solution z = might be stable before we guarantee that the evolution 
operator is a uniform contraction. In this case, however, we don't guarantee that 
stability of the trivial solutions persists under perturbations. Hence, we cannot guar- 
antee that the nonlinear perturbation coming from the Taylor remainder does not 
destroy the stability. We avoid tackling this case, since it would bring only further 
technicalities. Note the above a c synchronization is stable under small perturbations 



Example 3 Consider the Lorenz system presented in Sec. \3.2.1\ 
Then 

(-a -a a 
r-z —a- 1 -x 
y x —b — a 

noting that the trajectories lie within the absorbing domain Q given in Proposition 
[7] we have 



k\<y/r-i=, \y\<r r -^ = ! and|z-r|<r 



therefore, 



a t =r\ 1 + . + ^r—= - 1 

V s/o(b-\)) Vb^l 



For the standard parameters ( see Sec. \3.2.1\ we have (X c ~ 56 .2 . For the two coupled 
Lorenz, this provides the critical parameter for synchronization 

a > a c /2«28.1 



We have simulated the dynamics of Eq. ( |5.10 1 using the Lorenz system. For a 



27 > a c we observe that the complete synchronized state is stable. If the two Lorenz 
systems start at distinct initial condition as time evolves the difference vanishes 



exponentially fast, see Fig 5.1 



If we depict x\ x X2 the dynamics will lie on a diagonal subspace x — y. If the 
initial conditions start away from the diagonal x = y the evolution time series will 



then converge to it, see Fig 5.2 



5.4 Network Global Synchronization 



We turn to the stability problem in networks. Basically the same conclusion as be- 
fore holds: the network is synchronizable for strong enough coupling strengths. In 



5.4 Network Global Synchronization 
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a) b) 




5 10 15 20 0,2 0,4 0,6 0,8 1 



time time 

Fig. 5.1 Time evolution of norm ||*i(f) — *2(f)ll f° r distinct initial conditions, a) For a = 0.3 we 
observe an asynchronous behavior, b) for a = 21 above the critical coupling parameter the norm 
of the difference vanishes exponential fast as a function of times, just as predicted by the uniform 
contraction. 




Fig. 5.2 Behavior of the trajectories in the projection x\ x^. a) for the coupling parameter a = 3, 
the trajectories are out of sync, and spread around, b) for a = 21 the trajectories converge to the 
diagonal line x\ = X2. Trajectories in a neighborhood of the diagonal converge to it exponentially 
fast. 
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such a case we want to determine the critical coupling in relation to the network 
structure. A positive answer to these question is given by the following 

Theorem 7 Consider the diffusively coupled network model 

n 

Xi =f(xi) + a £ Ajj(xj -Xj), 
j=i 

on a connected network. Assume that the isolated system has a Lyapunov function 
satisfying Assumption^with an absorbing domain Q. Moreover, assume that for 
a given time s > all trajectories are in a neighborhood of the synchronization 



manifold lying on the absorbing domain £2, and consider the a c given by Eq. {5.13 ), 
and A2 the smallest nonzero eigenvalue of the Laplacian. Then, for any 

a c 

the global synchronization is uniformly asymptotically stable. Moreover, the tran- 
sient to the globally synchronized behavior is given the algebraic connectivity, that 
is, for any i and j 

\\xi{t)-*j{t)\\<Me- {ah - a ^ 



The above result relates the threshold coupling for synchronization in contribu- 
tions coming solely from dynamics a c , and network structure A2. Therefore, for 
a fixed node dynamics we can analyze how distinct network facilitates or inhibits 
global synchronization. To continue our discussion we need the following 

Definition 5 Let j3 (G) be the critical coupling parameter for the network G. We say 
that the network G is better synchronizable than H if for fixed node dynamics 

P(G)<P(H) 

Recalling the general bounds presented in Theorem[2]we conclude that the complete 
network is the most synchronizable network. Furthermore, the following general 
statement is also true 

- For a fixed network size, network with small diameter are better synchroniz- 
able. Hence, the ability of the network to synchronize depends on the overall 
connectedness of the graph. 



Recall the results presented in table 2.2 and let denote a c denote the critical 



coupling parameter, the dependence of a c in terms of the network size can be seen 
table El 

The difficulty to synchronize a complete network decreases with the network 
size, whereas to synchronize the cycle increases quadratically with the size. 

Now we present the proof of Theorem|7] We must show that the synchronization 
manifold ^# is locally attractive. In other words, whenever the nodes start close 
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Table 5.1 Leading order dependence of P on the network size for the networks in Fig. '. 



2.2 



Network /3 

Complete - 
n 



n 2 

ring — 

6 2 



Star 



together they tend to the same future dynamics, that is, \\xj(t) —Xj(t)\\ —> 0, for any 
i and j. For pedagogical purposes we split the proof into four main steps. 

Step 1: Expansion into the Laplacian Eigenmodes. Consider the equations of 
motion in the block form 

— = F{X)-a{L®I m )X 

Note that since L is symmetric, by Theorem [9] there exists an orthogonal matrix O 
such that 

L = OMO*, 

where M = diag(Ao, X\ , . . . , A„) is the eigenvalue matrix. Introducing 

Y = col(y u y 2 ,--- ,y„) 
we can write the above equation in terms of Laplacian eigenvectors 

X=(0®I m )Y, 

n 

1=1 

For sake of simplicity we calljj = s, and remember that now note that v\ = 1 
hence 

X = l®s + U, 

where 

n 
i=2 

In this way we split the contribution in the direction of the global synchronization 
and U, which accounts for the contribution of the transversal. Note that if U con- 
verges to zero then the system completely synchronize, that is X converges to 1 ®s 
which clearly implies that 
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The goal then is to obtain conditions so that U converges to zero. 



Step 2: Variational equations for the Transversal Modes. The equation of mo- 
tion in terms of the Laplacian modes decomposition reads 

^ =F{X)~a{L®I)X 
at 

1®-^ + — =F(l®s + U)-a(L®I) (l®s + U), 

We assume that U is small and perform a Taylor expansion about the synchroniza- 
tion manifold. 

F(l®s + U) = F{l®s)+DF(l®s)U + R(U), 

where R(U) is the Taylor remainder = 0(||t/|| 2 ). Using the Kronecker prod- 

uct properties [10] and the fact that LI = 0, together with 

ds 

/(gi — = F(l®s)=l®f(s) 
dt 



and likewise 
and we have 



DF(l®s)U = [I n ®Df(s)]U, 



^ = [DF(l®s)-a(L®l)]U+R(U), (5.14) 
dt 

Therefore, the first variational equation for the transversal modes reads 

= [/„ ® Df(s) -aL® I m ] U. 

The solution of the above equation has a representation in terms of the evolution 
operator 

U(t)=T(t,s)U(s) 

We want to obtain conditions for the trivial solution of the above to be uniformly 
asymptotically stable, that is, so that the evolution operator is a uniform contraction. 

Step 3: Stabilization of the Transversal Modes. Instead of analyzing the full set 
of equations, we can do much better by projecting the equation into the transversal 
modes v, . Applying v* ®I m on the right in the equation for U, it yields 

v* ®l m ^pvt ® = v* ®I m f £ Vi®Df{s) yi - aliVi ®y\ 

- dv- - 

i=2 al i=2 
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But since v, form an orthonormal basis we have 

v ) v i = 

where is 5;; the Kronecker delta. Hence, we obtain the equation for the coefficients 

d ^ = [Df{s)-al i l m ]y i 

All blocks have the same form which are different only by A,-, the z'th eigenvalue 
of L. We can write all the blocks in a parametric form 

d ^=K{t)u, (5.15) 

where 

K(t)=Df(s(t))-Kl m 

with K g K. Hence if K = aA, we have the equation for the z'th block. This is just 
the same type of equation we encounter before in the example of the two coupled 
oscillators, see Eq. ( |5.11| l. 

Now obtain conditions for the evolution operator of Eq. ( 5.15| l to possess a uni- 



form contraction. This is done applying the same arguments discussed in Eqs. |5.12| 



and 5.13 Therefore, the z'th block has a uniform contraction if aA, > a c . Now since 
the spectrum of the Laplacian is ordered, the condition for all blocks to be uniformly 
asymptotically stable is 

a c 

a > T - 
Aa 

which yields a critical coupling value in terms of a c and A2. 

Taking a larger than the critical value we have that all blocks have uniform con- 
tractions. Let Tj(t,s) be the evolution operator of the z'th block. Then 

\\yi(t)\\<\m,s)yi(s)\\ 

<\m, S )\\\\ yi ( S )i 

(5.16) 

by applying Theorem [5] we obtain 

l|y/(OII<^~'" M) l|y/(*)ll. 

where 77, = a A,- — a c . 

Step 4: Norm Estimates. Using the bounds for the blocks it is easy to obtain a 
bound for the norm of the evolution operator. Indeed, note that 



36 



5 Stability of Synchronized Solutions 



\u\\ 



1=2 



<Ll|Vi|IW2 

i=2 



where we have used Theorem 15 (see Ap.[A]i, therefore, 

\\U\\ 2 <t\\v,\\K,e-^'- s) \\y i (s)\\ 

1=2 

Now using that g~1;('~*) < e -( a h-"c)^ an( j applying Theorem|4]we obtain 

||T(M)l|2<M* -,,(f ~ ,) 

with tj = C1X2 — OL c for any f > s. 

By the principle of linearization Lemma[T] we conclude that the Taylor remainder 
does not affect the stability of the trivial solution, which correspond to the global 
synchronization. 

The claim about the transient is straightforward, indeed note that 

||ir(0-2®*(f)|| <Me-^'- s \\U{s)\\ 
implying that \\xi(t) -s(t)\\ < Ke^ f ^ and 

11^(0-^(011 < ||*i(0-*(0ll + ll*/(*)-*(0ll 

in virtue of the triangular triangular inequality, and we concluding the proof. □ 



Chapter 6 

Conclusions and Remarks 



/ have had my results for a long time: but I do not yet know how 
I am to arrive at them. 

— Gauss 



We have used stability results from the theory of nonautonomous differential equa- 
tion to establish conditions for stable global synchronization in networks of dif- 
fusively coupled dissipative dynamical systems. Our conditions split the stability 
condition solely in terms of the isolated dynamics and network eigenvalues. 

The condition associated with the dynamics is related to the norm of the Jacobian 
of the vector field. This reflects that fact that to obtain stable synchronization we 
need to damp all instabilities appearing in the variational equation. The network 
condition is given in terms of the graph algebraic connective - the smallest nonzero 
eigenvalue, which reflects how well connected is the graph. 

The dependence of synchronization on only two parameters is due to our hy- 
potheses: i) all isolated equations are the same, and if) the coupling between them 
is mutual and fully diffusive. These assumptions reduce the intrinsic difficulties of 
the problem and allow rigorous results. 

There are other approaches to tackling the stability of the global synchronization. 
Successful approaches are the construction of a Lyapunov function of the synchro- 
nization manifold, see for example Refs. Il46ll47ll48l , which takes a control view; 
and the theory of invariant manifolds |49l l50l taking a dynamical system view of 
synchronization. 

Our results have a deeper connection with the previous approach introduced by 
Pecora and Carrol ll26l . They used the theory of Lyapunov exponents, which allow 
the tackling of general coupling functions. The main drawback is that of obtaining 
results for the persistence of the global synchronization. This requires establishing 
results on the continuity of the Lyapunov exponent, which is rather subtle and intri- 
cate EQ.Q 

The approach introduced in these notes follows the steps of the Pecora and Carrol 
analysis, that is, the local stability analysis of the synchronization manifold, but uses 
various concepts in stability theory, to establish the persistence results for the global 
synchronization. 



1 Small perturbations can destabilize a system with negative Lyapunov exponents. To guarantee 
the persistence under perturbations Lyapunov regularity is required, see Ref. |51 1. 
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Appendix A 
Linear Algebra 



If only I had the theorems! Then I should find the proofs easily 
enough. 

- Bernard Riemann 

For this exposition we consider the field F where F = M. the field of real numbers or 
F = C the field of complex numbers. We shall closely follow the exposition of Ref. 
Il52l . Consider the set Mat(F, n) of all square matrices acting on F". We start with 
the following 

Definition 6 Let A G Mat( F, n ). The set 

(7(A) := {A £ C : det(A - Al) = 
is called the spectrum of A. 

The spectrum of A is constituted of all its eigenvalues. Note by the fundamental 
theorem of algebra the spectrum has at most n distinct points. 

Often, we want to obtain bounds on the localization of eigenvalues on the com- 
plex plane. A handy result is provided by the result 

Theorem 8 (Gershgorin) Let A e Mat(£,n), denote A = (Ay)? - =1 . Let D{a 7 8) 
denote the ball of radius 8 centered at a. For each 1 < i < n let 

Ri=i\Aij\, 

j=l 
m 

then every eigenvalue of A lies within at least one of the balls D{Au, R{). 
For a proof see Ref. [52] Sec. 10.6. 

If A 6 Mat(C,n) then we denote its conjugate transpose by A*. In case A is a 
real valued matrix A* denotes the transpose. A matrix is called hermitian if A = A* . 
The following definition is also fundamental 

Definition 7 Let A £ Mat(C,n) be a hermitian matrix. It is called positive-semidefinite 
(or sometimes nonnegative -definite) if 

x*Ax > 
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A. 1 Matrix space as a normed vector Space 
for any x G C" 

It follows that a matrix is nonnegative if all its eigenvalues are non negative. 
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A.l Matrix space as a normed vector Space 



Consider the vector space F" over the field F. A norm | ■ | onF" is a function 

|| -|| : F n -)• R satisfying 

1. positive definiteness : ||jc|| > for all* 6 F" and equality holds iff x = 

2. Absolute definiteness : \\yx\\ = \y\\\x\\ for all a: <G F n and ye F 

3. Triangle inequality : |*+j>| < ||xj| + ||3>| for all x,y e F" 

We call the pair (F n , \\ ■ ||) is called normed vector space. This normed vector 
space is also a metric space under the metric d : F" x F" — > R where d(x,y) = 
\\x—y\\. We say that d is the metric induced by the norm. In this metric, the norm 
defines a continuous map from F" to R, and the norm is a convex function of its 
argument. Normed vector spaces are central to the study of linear algebra. 

Once we introduce of norm on the vector space F", we can also view the 
Mat(F,n) as a normed spaces. This can be done by the induced matrix norm which 
is a natural extension of the notion of a vector norm to matrices. Given a vector 
norm || • || on F", we define the corresponding induced norm or operator norm on 
the space Mat(F,n) as: 

||A|| = sup{||A*j| :xeF n and = 1} 

It follows from the theory of functions on compact spaces that ||A|| always exists 
and it is called induced norm. Indeed, the induced norm defines defines a norm on 
Mat(F,n) satisfying the properties 1-3 and an additional property 

||AB|| < ||A||||B|| for all A, Be Mat(F,n) 

called sub-multiplicativity. A sub-multiplicative norm on Mat(F, n) is called matrix 
norm or operator norm. Note that even though we use the same notation ||A|| for the 
norm of A, this should not be confused with the vector norm. 

Example 4 Consider the norm of the maximum || • ||oo on M". Given R" 9 x = 
(jci, • • • ,x„), the norm is defined as ||x|| = max,- |x,-|. Given a matrix A = (A,- 7 -)" ; - =1 
then 

n 



max 

\<i<n 



Example 5 Consider the Euclidean norm \\ ■ \\2 on R". Using the notation of the 
last example, we have 

||A|| 2 = Vp(AM0 
where p max (A*A) is spectral radius A* A. 
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Recall that two norms || ||' and || ||" are said to be equivalent if 

fl||AH'< ||A||"<A||A||' 

for some positive numbers a,b and for all matrices A. It follows that in finite- 
dimensional normed vector spaces any two norms are equivalents. 



A.2 Representation Theory 



We review some fundamental results on matrix representations. A square matrix A 
is diagonalizable if and only if there exists a basis of F n consisting of eigenvectors 
of A. In other words, if the F" is spanned by the eigenvectors of A. If such a basis 
can be found, then P l AP is a diagonal matrix, where P is the eigenvector matrix, 
each column of P consists of an eigenvector. The diagonal entries of this matrix 
are the eigenvalues of A. One of the main goals in matrix analysis is to classify the 
diagonalizable matrices. 

In general diagonalization will depend on the properties of F such as whether 
F is a algebraically closed field. If F = C then almost every matrix is diagonaliz- 
able. In other words, the set B C Mat (C,«) of non diagonalizable matrices over C 
has Lebesgue measure zero. Moreover, the set diagonalizable matrices form a dense 
subset. Any non diagonalizable matrix, say Q G B can be approximated by a diago- 
nalizable matrix. Precisely, given e > there is a sequence {A,} of diagonalizable 
matrices such that \\Q — A;|| < £ for any i > no- 
Let us denote by * the conjugate transpose if F = C (clearly only transpose if 
F = M). We first focus on symmetric matrices A = A* and F = R. It turns out that 
it is always possible to diagonalize such matrices. 

Definition 8 A real square matrix A is orthogonally diagonalizable if there exists 
an orthogonal matrix P such that P* AP — T) is a diagonal matrix. 

Diagonalization of symmetric matrices is guaranteed by the following 

Theorem 9 Let A be a real symmetric matrix. Then there exists an orthogonal 
matrix P such that : 

1. P*AP = D is a diagonal matrix. 

2. D = diag (h\ , • • • , A„), where Xi are the eigenvalues of A.. 

3. The column vectors o/P are the eigenvectors of the eigenvalues of A. 



For a proof see Ref. 113311 Sec. 8.1. 



A.3 Kronecker Product 

A.3 Kronecker Product 
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We need several properties of the Kronecker Product to address the stability of the 
synchronized motion in networks. 

Definition 9 Let A G Mat ( F,mxn) and B e Mat( F,rxs). The Kronecker Product 
of the matrices A and B and defined as the matrix 



/AnB ••• Ai„B 



A(g)B: 



\A ml B 



AijiiiB 



The Kronecker product is sometimes called tensor product. Consider now the 
following examples on the 

Example 6 Consider the matrices 



A = 



a b 
c d 



and B = 



1 

2 3 



Then 



A<8>B : 



Now consider the vectors 



aB bB 
cB dB 



f a b 0\ 

2a 3a 2b 3b 
c Oa d 
\2c 3c 2d 3d J 



and u(t) 



x(t) 

y(t) 



Then 



V® u(0 



We review the basic results we need. 



fx(t)\ 

y{t) 

x(t) 



Theorem 10 Let A G Mat(F,m x nj one/ B e Mat(F, rx s) C e Mat(F,n x /?) anrf 
D € Mat(F,s x fj. r/ie« 

(A (g) B) (C (g) D) = AC (X) BD. 

The proof can be found in Ref. [52] pg. 408, see Proposition 2. Note that (A <E> 
B)(C<g)D) E Mat (F,mr x pt). A direct computation leads to the following result 



Theorem 11 Let A G Mat(F,m xnj andB G Mat(F,rx s), then 



(A(g)B)* = A* <g)B* 
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By applying Theorem[lO]we conclude that following 
Theorem 12 If A and B are nonsingular, then 

(A(K>B) _1 =A~ 1 <g)B- 1 . 

We following Theorem also plays a important role in the exposition 

Theorem 13 Let {A,}; =1 be the eigenvalues of A £ Mat(F,n) and {Hi} s i=1 be the 
eigenvalues of B £ Mat(F,n). Then A®B /zas rs eigenvalues 

h Mi , Ai/x,, , • • • , A2ju s , •», ^-Kv ■ 

The proof can be found in Ref. [52] pg. 412. A direct consequence of this result 
is the following 

Theorem 14 Let A and B be positive semi-definite matrices. Then A ® B is also 
positive semi-definite. 

Our last result concerns the norms of the Kronecker products 

Theorem 15 Let \\ ■ |L be p-norm. Consider v £ W, andx £ M?,fort,s £ N. Then 

h®x\\ P = h\\ P \M\p 



Appendix B 

Ordinary Differential Equations 



Let D be an open connected subset of W, m > 1, and let G : D — > WL m be an au- 
tonomous vector field. Consider the problem of finding solutions for the vector dif- 
ferential equation 

x = G(x) (B.l) 

with the initial condition x(0) =xq. A positive answer to this problem is given by 
the following 

Theorem 16 (Picard-Lindelof) Assume that the vector field G Lipschitz continu- 
ous in a neighborhood of xq. Precisely, assume that given xo 6 U C D there is a 
constant Kfj such that 

||G(x)-G(y)|| <Ku\\x-y\\ 



for all x, y € U. Then there exists a unique local solution x(f) for Eq. (B.l \ satisfy- 
ing x(0) =x . 

Note that the solution is local, in the sense that there is small K > such that the 
function* : [— K, fc] — > D is a solution of the problem wifhjc(O) =xq. The question is: 
How long does such solution exist for? We are interested in the long term behavior 
of the solutions, so we wish to know under what conditions the solutions exists 
forward in time. A positive answer is given by extension theorems: 

Theorem 17 (Extension) Let ^ be a compact subset of the open set D. Consider 



Eq. (B.l \ and let G be differentiate. Let Xo G c € and suppose that every solution 
x : [0, t] — > D with x(0) = x/o lies entirely in C. Then this solution is defined for all 
(forward) time t > 0. 

The proofs of the above theorems can be founds in Refs. Il54ll55l . 
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44 B Ordinary Differential Equations 

B.l Linear Differential Equations 



The evolution operator also determines the behavior of the non homogeneous equa- 
tion 

Theorem 18 Let A : M — > Mat(R,n) and g : M. — > W be continuous function. Con- 
sider the perturbed equation 

y = Ay + g(r) 

The solution of the perturbed equation corresponding to the initial condition x(fo) = 
xo is given by 

y(?)=T(f,fo)x + f'j(t,s)g(s)ds 

where T(f,fo) is the evolution operator of the corresponding homogeneous system. 

The following inequality is central to obtain various estimates 

Lemma 2 (Gronwall) Consider U C R+ and let u:U — > Ml Z?e continuous and non- 
negative function. Suppose there exist C > and and K > such that 

u(t) <C+ f Ku(s)ds (B.2) 
Jo 

for all t GU, then 

u(t) <Ce Kt . 

The proof of these results can be found in Ref. Il54l . 
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